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ABSTRACT 

The  quantum  Boltzmann  equation  method  is  demonstrated  by  numerically  predicting  the  time-dependent  so¬ 
lutions  of  the  velocity  and  magnetic  fields  governed  by  nonintegrable  magnetohydrodynamic  equations  in  one 
spatial  dimension.  The  method  allows  arbitrary  tuning  of  the  value  of  the  viscosity  and  resistivity  transport 
coefficients  without  compromising  numerical  integrity  even  near  the  zero  dissipation  and  turbulent  regime  where 
shock  front  discontinuities  emerge. 

Keywords:  quantum  computing,  quantum  Boltzmann  equation,  magnetohydrodynamic  turbulence 

1.  INTRODUCTION 

Parallel  quantum  algorithms1  have  been  developed  to  numerically  predict  time-dependent  field  solutions  of  the 
diffusion  equation,2  the  nonlinear  Burgers  equation,3,4  the  one  dimensional  equations  for  magnetohydrodynamic 
turbulence,5  the  Korteweg  de  Vries  equation,6  the  nonlinear  Schroedinger  equation,6,7  and  the  Manakov 
equations.8,9  The  first  parallel  quantum  algorithms  for  the  diffusion  and  Burgers  equations  have  with  success 
been  experimentally  tested  on  quantum  information  processing  prototypes  of  parallel  quantum  computers  (or 
type-II  quantum  computers10)  using  spatial  nuclear  magnetic  resonance  spectroscopy  on  a  linear  array  (in  space 
and  in  reciprocal  space)  of  segmented  ensembles  of  two-qubit  labeled  chloroform  molecules.11-13 

Here,  we  generalize  our  parallel  quantum  algorithm  for  modeling  magnetohydrodynamic  (MHD)  turbulence. 
This  algorithm  is  based  upon  the  quantum  algorithms  for  the  diffusion2  and  Burgers  equations3,4;  it  can  be 
implemented  on  a  parallel  quantum  computer  with  pairs  of  labeled  chloroform  molecules.  The  intent  of  this  pa¬ 
per  is  solely  to  demonstrate  the  quantum  Boltzmann  equation  method,  so  all  the  details  about  the  microscopic 
quantum  mechanical  implementation  of  the  method  have  been  replaced  with  a  friendlier  mesoscopic  treatment. 
The  quantum  mechanical  foundation  of  the  method  lends  it  several  advantageous  features  as  a  tool  for  compu¬ 
tational  physics.  First,  there  is  inherent  numerical  stability.  Second,  there  is  a  corresponding  entropy  function 
description  of  the  model  dynamics  and  an  entropy  theorem,  although  for  economy  and  brevity  the  details  of  the 
entropy  function  are  not  presented  here.  The  Boltzmann  collision  function  can  be  computed  directly  without  the 
need  of  root  finding  as  is  required  for  entropic  lattice  Boltzmann  equation  models  of  turbulence.14-16  Although 
the  quantum  algorithm  presented  here  was  originally  developed  for  a  parallel  quantum  computing  architecture, 
the  quantum  Boltzmann  equation  is  a  practical  algorithmic  method  for  presentday  classical  computers. 

Turbulence  is  plagued  by  spatiotemporal  intermittency  involving  coherent  structures  -  structures  that  are 
at  odds  with  the  simple  scale-similarity  arguments.  Insight  into  these  coherent  structures  comes  by  examining 
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simplified  models  and  the  simplest  model  exhibiting  Alfvenization  (the  exchange  of  fluid  and  magnetic  energies) 
are  the  ID  equations17,18 


dv  8  ( B2\  d2v 

+  Vdx  =  ~dx\ir)  +  f*!h? 


dB  d  .  d2B 

ir  +  ai(«s>  =  ’a?  (2) 

for  the  velocity  field  v  =  v(x)x  and  the  magnetic  field  B  =  B(x)z.  The  transport  coefficient  p  is  the  fluid 
viscosity  and  77  is  the  plasma  resistivity.  These  ID  equations  are  derived  from  the  full  MHD  equations  when  the 
fluid  density  length  scales  are  much  longer  than  those  of  the  B  field.17, 18  To  leading  order,  this  results  in  a 
constant  density  so  that  (1)  forms  a  closed  set  of  equations  for  v  and  B.  Adding  a  ^-component  of  the  magnetic 
field  gives  a  ID  model  for  solar  flares.19  In  Elsasser  variables  z±  =  v  ±  B,  (1)  can  be  written  in  symmetric  form 


dz±  dz± 
dt  2±  dx 


•±5(9-/*) 


d2(z+  -  zJ) 


If  the  transport  coefficients  are  equal,  77  =  p,  then  (3)  reduces  to  two  uncoupled  nonlinear  Burgers  equations 
for  z±.  [In  the  limit  B  =  0,  (1)  reduces  to  Burgers  equation  for  v].  The  Burgers  equation20  is  a  well-known 
paradigm  for  Navier-Stokes  turbulence  and  has  also  been  studied  extensively  for  many  decades  as  a  simplified 
model  for  boundary  layer  behavior,  shock  wave  formation,  mass  transport,  self-organized  criticality  and  growing 
interfaces.  It  is  also  a  test-bed  for  numerical  methods  since  a  general  analytic  solution  exists.21,22  In  Burgers 
turbulence,  regions  where  <  0  steepen  into  shock  singularities,  while  regions  where  >  0  become  smoother. 
However,  with  the  inclusion  of  the  B  field,  there  is  now  a  magnetic  back-pressure  in  (1)  that  is  absent  from 
Burgers  equation.  In  particular,  the  B  field  will  concentrate  in  regions  of  the  velocity  shock,  softening  the  shock 
front.  For  77  7^  77,  (1)  is  non-integrable. 

To  model  (1)  we  use  two  quantum  Boltzmann  equations  in  a  two  step  algorithm  for  an  ordered  one  dimensional 
array  of  L  number  of  nodes,  enumerated  by  a  “spatial”  integer  coordinate  x  =  1,2, ...,  L.  At  each  node  are  four 
probability  values  0  <  pa  <  1,  for  a  =  1,2, 3, 4  and  they  are  called  occupation  probabilities.  The  first  step  of  the 
algorithm  encodes  the  B[x)  field  using  the  probability  fields  p3(x)  and  pi(x): 


Ps(x)  = 


Pi(x)  = 


The  second  step  uses  the  following  quantum  Boltzmann  equation  to  update  these  two  probability  fields: 

p'3{x-l)  =  p3  {x)  -  [p3  (x)  -  pi  (as)]  sin2  /? 

Pi(x  +  1)  =  Pi{x)+  \p3(x)  -p4(x)]sin2^. 

The  third  step  assigns  the  new  values  to  the  B(x)  field: 

B(x)  =p'3(x)  +p'4(x). 

The  fourth  step  encodes  the  v(x)  and  B(x)  fields  using  all  four  probability  fields: 


Pi,2(as)  =  -[i  +v(x)  ±B{x)\  =p3A(x). 

(Only  on  the  first  pass  through  the  algorithm,  to  initialize  it  in  local  equilibrium,  we  use  the  encoding: 

Pi  (x)  —  di  [cot  9, 1  4-  v(x)  +  B(x)] 

P2(x)  =  d2{cot9, 1  +  v(x)  -  B(x)\ 

Ps{x)  =  di[cot0, 1  +  v(x)  -f  B(x)] 

Pi{x)  =  ^[cotfl,  1  +  v(x)  —  f?(x)], 
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where  the  local  equilibrium  functions 

di(a,p)  =  |  +  ^  [\/l  +  a2-  Vl +«2(p- l)2]  (13) 

&(<*,?)  =  2  ~  ^  [\/l  +  a2-  -/l +a2(p- 1)2J  (14) 

are  the  zeros  of  the  quantum  collision  function.)  The  fifth  step  uses  the  following  quantum  Boltzmann  equation 
to  update  all  the  probability  fields: 

p[{x  -  1)  =  Pi{x)  -  \pi{x)-p2{x)}s\ne  +  y/pi{x){l-pi(x)\p2{x)[l-p2{x)\  sin 20  (15) 

p'2(x  +  l)  =  p2{x)  +  [pi(z) -p2{x)\smO  -  -\/pi (^)[1  -  Pi(x))P2{x)[1  ~P2(x)}  sin 20  (16) 

p'3{x-  1)  =  p3(x)  -  [p3(a:)  -p4(a:)]sin0+  v/p3(a:)[l  -P3(a:)b4(a:)[l  -P4(a:)]  sin  20  (17) 

p'4(x+l)  =  Pi(x)  +  \p3(x)  -  p<(x)\  sm9  -  (18) 

The  final  step  uses  an  inverse  of  (8)  to  compute  the  updated  values  of  the  v(x)  and  B(x)  fields  from  the  new 
occupations  probabilities  by  the  following  assignments: 

v'(z)  =  )^\pi{x)  +  p2{x)  +  p3{x)  +  pi{x)\  -  l  (19) 

B'  (x)  -  i  [pi  (x)  +  p2  {x)  -  p3  (x)  -  p4  (x)] .  (20) 


At  this  point,  one  iteration  of  the  algorithm  is  completed  and  time  is  advanced  by  one  unit  At.  In  this  way,  a 
time  history  of  the  v  and  B  fields  can  be  generated  by  looping  through  the  enumerated  steps. 

2.  EFFECTIVE  FIELD  THEORY 

After  some  algebra  for  the  special  case  where  /3  =  0  and  0  =  f ,  and  with  the  encoding  (8)  and  the  asignment 

(1)  formulae,  it  follows  that  the  mesoscopic  transport  equation  (1)  can  be  written  directly  as  a  finite  difference 
equation  in  the  v(x)  and  B(x)  fields: 

v'(x)  =  ~  [u(x  -  1)  4-  v(x  +  1)]  +  ^  [v(x  -  l)2  -  v(x  +  l)2]  +  i  [B(x  -  l)z  -  B(x  +  l)2]  (21) 

B'{x)  =  i  [B(x  +  1)  +  B(x  -  1)]  +  i  [v(x  -  \)B(x  -  1)  -  v(x  +  \)B{x  +  1)] .  (22) 

A  more  general  set  of  finite  difference  equations  than  (2)  can  be  predicted  by  the  quantum  Boltzmann  equation 
method  for  arbitrary  9,  but  are  too  lengthy  to  include  here  *. 

The  macroscopic  scale  effective  field  theory  for  v(x,t)  and  B(x,  t)  follows  from  taking  the  continuum  limit  of 

(2)  with  the  ordering  Ax  ~  (D(e)  and  At  ~  0(e2).  We  obtain  (in  lattice  units  Ax  =  1  =  At) 


dv  dv 

1  d2y 

—  dB  ». 

Ht+Vdx 

2  dx2 

-  B  ^  +  0(«») 

(23) 

dB  dB 

1  d2B 

_ dv  0. 

dt  +Vdx 

2  dx2 

(24) 

This  is  just  the  ID  MHD  (1)  in  the  special  case  where  the  kinematic  shear  viscosity  and  resistivity  transport 
coefficients  are  equal,  p  =  r)  =  If  we  make  an  Elsasser  variables  substitution,  z±  =  v±B,  then  (23)  and  (24) 
simplify  into  two  uncoupled  nonlinear  Burgers  equations  for  z±,  see  (3).  We  now  proceed  from  the  integrable 
ID  MHD  with  p  —  r)%o  the  non-integrable  case  of  unequal  transport  coefficients,  p  ^  77. 

‘The  unitary  collision  operators  are  implemented  with  perfect  control  and  the  occupation  probabilities  are  measured 
exactly.  These  assumptions  hold  only  for  an  ideal  experimental  setup.  We  address  actual  quantum  control  and  Von 
Neuman  projective  shot  noise  in  our  papers  on  experimental  parallel  quantum  computing  with  NMR.11"13 
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After  some  algebra  for  the  special  case  where  ft  =  f  and  0  =  f ,  and  with  the  encoding  formulae  (4)  and  the 
mesoscopic  transport  equation  (1)  and  the  asignment  formula(7),  and  then  also  with  the  encoding  formulae  (8) 
and  the  quantum  Boltzmann  equation  (1)  and  the  asignment  formulae  (1),  it  follows  that  the  dynamics  of  the 
v(x)  and  B(x)  fields  can  be  written  directly  as  a  finite  difference  equation: 

v'(x)  —  i  [v(x  -  1)  +  v(x  +  1)]  +  ^  [v(x  -  l)2  -  v(x  +  l)2]  +  (25) 

4  [B(x  -  2)  -  B(x  +  2)]  [. B(x  -  1)  +  2B(x)  +  B{ x  +  1)] 

B'( x)  =  i  [B{x  -  2)  +  2B{x)  +  B(x  +  2)]  +  (26) 

i  [v(x  —  1  )B(x  —  2)  +  [v(x  —  1)  —  v(x  +  l)]f?(a;)  —  v(x  +  1  )B(x  +  2)] . 

Again  proceeding  to  the  continuum  limit  by  Taylor  expanding  (2),  we  obtain  our  desired  ID  MHD  equations  with 
unequal  transport  coefficients:  viscosity  p  =  ^£7  and  resisitivity  tj  =  In  the  general  case  where  ~  <  0  <  | 
and  f  <  P  <  f,  we  obtain  the  following  effective  field  theory: 


dv 

dv 

d2v 

dt  V  dx 

=  ^ 

dB 

dB 

d2B 

dt  +V  dx 

Bfx  +^3) 


with  adjustable  transport  coefficients  for  the  viscosity  and  resistivity  * : 


fi  =  cot2  0- 


rj  =  (cot2  9  +  cot2  (3) 


3.  NUMERICAL  TESTS 

The  ID  MHD  equations  are  solved  using  the  quantum  Boltzmann  equation  method  described  above  for  initial 
profiles  v(x,  0)  =  cos  (^)  and  B(x,  0)  =  O.lon  ai  =  256Ax  grid.  The  evolution  of  the  fields  are  shown  in  Fig.  1 
(low  dissipation  regime)  and  Fig.  2  (high  dissipation  regime).  Regions  of  <  0  steepen  towards  shocks  while 
regions  of  >  0  smooth.  At  first  the  B  field  rapidly  increases  by  a  factor  of  10  and  concentrates  where  velocity 
shocks  form.  Asymptotically  in  time  (but  when  the  nonlinear  terms  are  important  and  before  the  profile  linearly 
diffuses  away)  the  velocity  profile  bifurcates  and  eventually  becomes  linear  in  x  except  at  the  shocks.  There 
is  a  surprising  cancellation  of  errors.  If  you  look  at  the  profiles  of  the  occupation  probabilities  in  Fig.  1  (low 
dissipation  regime),  they  separately  each  have  Gibbs  oscillations.  However,  when  the  two  occupation  probability 
fields  are  added  together  forming  the  v(x)  field  or  subtracted  forming  the  B(x)  field,  the  oscillations  substantially 
cancel  out  because  the  oscillations  are  out  of  phase.  The  Gibbs  oscillations  subside  altogether  with  higher  grid 
resolution  since  its  source  is  an  under  resolved  shock  front. 

The  B  field  asymptotically  tends  to  a  constant  with  spatial  jumps  around  the  velocity  shock  fronts.  The  linear 
sloped  v  field  decays  in  time  as  does  the  constant  B  field.  Indeed,  we  can  check  our  model  in  this  asymptotic 
region.  Assuming  v(x,t)  — >  a(t)x  +  b(t)  and  B(x,t)  — >  B(t),  except  at  isolated  locations  corresponding  to 
shock  fronts,  (1)  yield  constraint  equations  for  a(t),  c(t),  and  B(t): 

da  o  db  dB 

~77  +  u  =0  —  +  af>  =  0  “-7- — |-  aB  =  0.  (30) 

dt  dt  dt 

The  analytical  solution  of  (30)  yield  asymptotic  v  and  B  field  profiles 


v(x,t)  = 


X  +  k2 


B(t)  = 


+  (29)  based  on  previous  results  from  modeling  the  Burgers  equation,  but  their  exact  form  has  not  yet  been  confirmed. 
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Figure  1.  Simulation  results  in  the  low  dissipation  regime  /?  =  0  and  6  =  1.5  radians.  The  blue  curve  is  the  velocity  field. 
The  red  curve  is  the  magnetic  field.  The  purple  curves  are  the  first  and  second  occupation  probability  fields  (pi  and  P2). 
The  green  curves  are  the  third  and  fourth  occupation  probability  fields  ( P3  and  pi).  There  are  four  black  curves;  these 
are  analytical  predictions,  based  on  the  computed  macroscopic  field  profiles,  of  the  equilibrium  occupation  probabilities 
and  they  are  in  excellent  agreement  with  the  numerical  predictions.  The  mesoscopic  fields  remain  in  local  equilibrium. 
The  simulation  was  carried  out  on  a  lattice  with  256  nodes.  Snapshots  are  shown  at  every  512  time  steps.  In  regions  of 
<  0,  shocks  tend  to  form  in  the  u-field  with  localized  enhancement  of  the  B- field.  Asymptotically,  v  attains  a  constant 
slope  with  shock  discontinuities  while  the  B  field  is  constant  with  discontinuities  at  the  leading  edges  of  the  front.  The 
extent  of  the  Gibbs  oscillations  at  the  shock  front  is  greater  at  the  mesoscopic  than  at  the  macroscopic  scale  because  of 
a  fortuitous  cancellation  of  errors. 
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time 


Figure  3.  The  time  evolution  of  the  B  field  at  location  x  —  192  determined  from  the  quantum  algorithm  (dotted  curve). 
Two  cases  are  shown.  For  the  6  =  \  case,  after  the  shock  formation  around  t  ~  200Af .  For  the  8  —  1.5  radians  case, 
shock  formation  occurs  latter  in  time  at  around  t  ~  2500Af.  The  analytical  asymptotic  solutions  (red  curves)  shows  the 
B  field  decays  like  t~x  as  expected,  see  (31).  There  is  excellent  agreement  between  theory  and  the  numerical  solution  over 
the  whole  dissipative  regime  from  high  to  low  resistivity. 

with  constants  of  integration  k\ ,  k?  and  k3.  From  the  quantum  simulations  shown  in  Figs.  1  and  2,  we  plotted 
in  Fig.  3  the  B  field  at  x  =  192  as  a  time  series  (black  dots).  The  solid  curve  is  the  asymptotic  decay  (31) 
derived  from  (1).  The  decay  is  fast  for  the  high  resistivity  case  and  it  is  slow  for  the  low  resistivity  case.  There 
is  excellent  agreement  after  the  shock  has  fully  formed  with  the  analytical  asymptotic  solutions  (red  curves) 
validating  our  model  for  ID  MHD  with  unequal  and  arbitrary  transport  coefficients. 

4.  FINAL  REMARKS 

The  numerical  simulation  results  plotted  in  Figs.  1  and  2  manifest  a  characteristic  property  of  a  type-II  quantum 
algorithm.  The  expected  value  of  the  occupancy  of  the  ground  state  and  excited  state  of  a  microscopic  qubit 
contained  within  a  quantum  mechanical  node  of  the  lattice  (in  the  ID  MHD  model  presented  here  there  would 
be  four  qubits  per  node)  are  shifted.  The  expectation  value  of  the  excited  state  is  plotted  in  the  purple  and  green 
curves  shown  in  in  Figs.  1  and  2.  Regardless  of  the  dissipation  regime  (high  or  low  viscosity  and  resistivity), 
there  is  a  gap  in  value  of  the  occupation  probabilites;  hence,  p°i  ^  p"]'  (two  purple  curves)  and  p*3  /  p°£  (two 
green  curves).  The  physical  cause  of  this  gap  is  the  following.  The  microscopic  dynamics  is  interpreted  as  the 
motion  of  quantum  particles  moving  left  and  right  through  a  chain  of  quantum  processors.  If  the  likelihood 
of  a  quantum  particle  exiting  left  or  right  from  a  quantum  processor  is  equal,  then  the  excited  state  energies 
encoding  the  left-  and  right-going  probabilities  would  be  degenerate.  The  equilibrium  probabilities  must  overlap: 
VT  —  P°2  arK-l  P°3  —  PT-  This  degeneracy  occurs  in  pairs  because  there  are  four  occupation  probabilities  and 
only  two  macroscopic  field  quantities.  The  macroscopic  effective  field  theory  would  be  strictly  diffusive  as  any 
quantum  particle  would  move  up  and  down  the  one  dimensional  chain  of  quantum  processors  in  an  unbiased 
random  walk  fashion.  However,  if  the  likelihood  of  a  quantum  particle  exiting  a  quantum  processor  to  the  left 
does  not  equal  the  likelihood  it  will  exit  to  the  right,  then  the  degeneracy  in  the  distribution  functions  would 
be  lifted  and  an  energy  gap  would  appear  as  demonstrated  Figs.  1  and  2.  In  this  case,  the  macroscopic  effective 
field  theory  would  not  be  strictly  diffusive;  there  would  be  an  overall  net  advection  of  quantum  particles  in  one 
direction  (the  symmetry  of  the  lattice  is  broken  and  this  causes  an  energy  gap  in  excited  state  energy  levels). 
This  net  advection  gives  rise  to  the  nonlinear  terms  in  (2). 

The  quantum  Boltzmann  equation  is  a  new  modeling  tool  for  numerically  predicting  the  time-dependent 
behavior  of  magnetohydrodynamic  turbulence  and  shock  formation.  Further  tests  of  the  algorithm  will  be 
presented  in  a  subsequent  paper.  An  open  problem  is  the  extending  the  method  to  handle  magnetohydrodynamic 
turbulence  in  two  and  three  spatial  dimensions. 
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